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Abstract. If the diffuse extragalactic gamma ray emission traces the large scale 
structures of the universe, peculiar anisotropy patterns are expected in the gamma 
ray sky. In particular, because of the cutoff distance introduced by the absorption of 
0.1-10 TcV photons on the infrared/optical background, prominent correlations with 
the local structures within a range of few hundreds Mpc should be present. We provide 
detailed predictions of the signal based on the PSCz map of the local universe. We also 
use mock N-body catalogues complemented with the halo model of structures to study 
some statistical features of the expected signatures. The results are largely independent 
from cosmological details, and depend mostly on the index of correlation (or bias) of 
the sources with respect to the large scale distribution of galaxies. For instance, the 
predicted signal in the case of a quadratic correlation (as it may happen for a dark 
matter annihilation contribution to the diffuse gamma flux) differs substantially from 
a linear correlation case, providing a complementary tool to unveil the nature of the 
sources of the diffuse gamma ray emission. The chances of the present and future space 
and ground based observatories to measure these features are discussed. 
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1. Introduction 

Gamma ray astronomy is a nourishing field in astroparticle physics. Results have 
rapidly accumulated in the last decade or two: after the break-through results of the 
EGRET satellite, a whole series of Earth-based observatories have developed, both using 
the imaging air Cherenkov telescopes (ACT like WHIPPLE, HEGRA, CANGAROO, 
HESS, and MAGIC) or surveying the sky via extensive air showers (EAS detectors, like 
TIBET, ARGO, and MILAGRO). New projects using all these techniques (the ACTs 
VERITAS, HESS II, MAGIC II, MACE, the EAS detector HAWC, the satellites AGILE 
and GLAST) are on the way (for a recent review of the field, see PQ). 

At TeV energies, a few dozen sources have been detected mainly by air Cherenkov 
experiments, most of which are high energy counterparts of MeV-GeV sources in the 
EGRET catalog. The 0.1-10 TeV range represents one of the "last" photonic windows 
yet to be explored at large distances. Starting from an energy of about 100 GeV (which 
we shall refer to as very-high energy, VHE), the absorption of high energy photons onto 
the extragalactic background light (EBL) via pair-production introduces an energy- 
loss horizon of the order of a few hundreds Mpc or smaller, well below the size of the 
observable Universe. This horizon drops down to ~10 kpc at PeV energies, basically 
precluding deep space astronomy with photons above about 10 TeV. Far from being 
only a limitation, this phenomenon also allows the use of 7-ray astronomy to probe the 
EBL, which is otherwise difficult to study directly. 

Besides single sources, wide field of view instruments and satellite-based 
observations are sensitive to diffuse 7-ray emissions. A particularly interesting 
emission is the extragalactic diffuse 7-ray background (in the following, cosmic gamma 
background, or CGB). The CGB is a superposition of all unresolved sources emitting 7- 
rays in the Universe and provides an interesting signature of energetic phenomena over 
cosmological time-scales. While a clear detection of this background has been reported 
by the EGRET mission [2], its origin is still uncertain, despite the fact that many 
models have been proposed. The most likely contribution is the one from unresolved 
blazars, i.e. beamed population of active galactic nuclei [3], with (probably sub-leading) 
components from ordinary galaxies [I], clusters of galaxies [5], and gamma ray bursts [6]. 
However, exotic possibilities like dark matter annihilation have been proposed, that are 
compatible with existing data and constraints [3, [HI El ID]. It is extremely difficult 
to test such models as long as the only observable is the energy spectrum. Recently, 
it was proposed to use the peculiar small-scale anisotropy encoded in the MeV-GeV 
gamma sky to probe dark matter [H] or astrophysical (T2J [13] contributions to the 
CGB. In this work, we further study this topic, with particular emphasis on the large 
scale anisotropy in the energy range 0.1-10 TeV. The lower part of this range will be 
probed by the GLAST telescope [Ml [15], while the energy window above the TeV is in 
principle accessible to EAS detectors. Different candidates to explain the CGB predict 
distinctive large scale features, even when similar energy spectra are expected. This is 
a consequence of the combined effect of a cutoff distance after which VHE 7 can travel 
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undamped to us, and of the anisotropic distribution of matter in the local universe (i.e., 
within a few hundred Mpc from us). A similar anisotropy pattern in the the ultra-high 
energy cosmic ray sky was recently analyzed by some of the authors [16] (see also [17]). 
The main goal of this paper is to characterize various features of the CGB which may 
be used for diagnostics at VHE. 

We do not attempt here to derive the properties of the CGB from specific 
astrophysical models, such as blazars (see e.g. [18]) or large-scale structure shocks 
(see [T9[ 1201 121]). Instead, we shall consider a phenomenological and parametric 
approach, analyzing in detail the angular patterns vs. energy for two representative 
diffuse background models: in the first case, we assume that underlying gamma emitters 
correlate linearly with overdensities in the large scale structure density field. In the 
second case, we assume a quadratic correlation. Loosely, one may consider the first 
case as representing unbiased, unresolved astrophysical sources; the second one may be 
indicative of a strongly biased astrophysical population of sources, or eventually of dark 
matter annihilation emission^ provided it is not dominated by substructures (see [22] for 
a discussion of this point). Of course, the yet unknown sources of the CGB may have a 
non-trivial bias and evolution with redshift. To a great extent, additional information on 
this topic will be provided by the study of the numerous point-like sources that GLAST 
is expected to resolve. At present, our analysis should be considered as a parametric 
approach to evaluate the sensitivity to "effective" matter tracing properties and bias of 
gamma ray sources. 

The paper is structured as follows. The parametrization we use for the CGB, 
based on EGRET data, is summarized in Sec. [21 In Sec. [3] we introduce our treatment 
of the cosmic large scale structure: a mock catalogue derived from a dark matter N- 
body simulation and the 3D power spectrum derived from the Halo Model of non- 
linear clustering J23[ [231 I2S1 [36]. With these tools, we shall argue that in the VHE 
range, especially at the largest angular scales, the predictions mostly depend on the 
large scale structure in the local neighborhood of the Universe. We shall then use the 
PSCz astronomical catalogue as tracer of the real structures in the nearby universe, 
thus producing maps of the VHE gamma sky (Sec. [3J. In Section [5] we analyse the 
perspectives for detection of these features using the forthcoming satellite GLAST, 
and for EAS observatories like MILAGRO. In Sec. [6] we summarize our findings, and 
conclude. In Appendix A we provide some details on the parametrization of the EBL, 
and on the method we use to propagate gamma rays and calculate attenuation effects. 



Appendix B reviews some statistical properties of a discrete poisson process on the 



sphere, relevant for our estimates of the errors of the multipole maps. 

| Of course, dark matter emission would also have peculiar features in the energy spectrum (e.g. 
departure from a power-law shape), associated with the particle physics details. They may constitute 
an important diagnostic tool for detection. 
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2. The cosmic gamma background 

Experimentally, the CGB is the most difficult component of the diffuse emission to 
study. Indeed, it is not correct to assume that the isotropic component after extracting 
point-like sources plus the galactic diffuse emission is entirely extragalactic: even in 
the pole direction, the CGB is comparable to the Galactic contribution. The deduced 
CGB thus depends on the adopted model of the Galactic foreground. The analysis 
undertaken to derive the spectrum of the CGB based on EGRET data provided the 
intensity spectrum [2] 

/ E7 \ -2.10±0.03 

valid from E ~ 10 MeV to E ~ 100 GeV, where k = (7.32 ± 0.34) x 
10 -6 cm -2 s -1 sr -1 GeV - . Interestingly, it shows a spectral index remarkably close to the 
average one of 7-ray blazars detected by EGRET, —2.1 ±0.3 [26]. EGRET experimental 
points and the best fit curve are shown in Fig. [TJ 

It is worthwhile to comment that the foreground subtraction remains a delicate 
issue, as can be appreciated by the reanalysis of the data performed in [27] , based on 
a revised model for the galactic propagation of cosmic rays. The deduced extragalactic 
spectrum is significantly lowered with respect to Eq. ([I]) at intermediate energies, while 
closer to the original result of Eq. (TjQ) at the lowest and highest energy points. Since 
the removal of the isotropic galactic emission from the diffuse gamma background is an 
open problem, and we still lack a complete understanding of the sources of the CGB, 
we shall simply base our following analysis on the extrapolation of the spectrum of 
Eq. (TjQ) by one to two orders of magnitude. Actually, the anisotropy pattern in the 
CGB sky itself may help in the foreground removal. Recently, it was proposed to use 
the cosmological Compton-Getting effect (whose dipole direction and amplitude are 
basically energy- independent) to discriminate the truly extragalactic fraction from the 
galactic foreground [28] . 

When extrapolating the EGRET flux to higher energies attenuation effects must 
be taken into account. It is well known that the propagation of photons in the EBL is 
a crucial issue for gamma astronomy in the VHE range. Absorption of VHE photons 
through pair-production with CMB, infrared or optical photons distorts an initial source 
spectrum, in particular by steepening its high energy tail. Astronomy of gamma ray 
emitters like blazars and active galaxies requires then an accurate modelling of the 
photon propagation and of the background frequency distribution. In Appendix A we 
provide a detailed description of the model of the EBL used, and of the treatment of 
absorption effects. For an overview of the knowledge of the EBL we refer the reader to 
the review [29]. The expression for the flux expected at Earth in a generic cosmology is 
given by (see e.g. [30l IB]) 



00 



HE n)oc I dz (2 ) 
I(E„n)ctJ Q dz H(z)(l + z)* ' (2) 
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Figure 1. EGRET spectrum from [5] and extrapolation up to 10 TeV. The dashed 
line shows the expected effect of the pair-production attenuation. 



where we assume an universal spectrum for the source, g(E), E y is the energy we observe 
today, p(z, n) is the matter density in the direction n at the redshift z, where the sources 
are assumed to be distributed proportionally to p a . The Hubble parameter is related to 
the present Hubble expansion rate H through the matter and the cosmological constant 



energy density as H(z) = H JQ M (1 + z) 3 + Q\. The quantity t(E^,z) is the optical 
depth of photons to absorptions via pair production on the EBL (see Appendix A). For 
most of the following considerations the normalization in the spectrum is irrelevant. 
Where ever needed (e.g. to estimate the statistics which can be collected by a given 
experiment) we shall fix the normalization of Eq. (jSj) so that it matches the EGRET fit 
of Eq. p at 10 GeV. 

For most of what follows it is important to realize that, by looking at VHE, most 
of the dependence on the cosmology, the source evolution, etc. in Eq. (J2J) cancel 
out, independently of the index of correlation with density, because of the cutoff at 
z 1 existing for VHE gammas. This important property dramatically reduces 
the model-dependence of the following considerations. It is also worth commenting 
that the 7-rays of the CGB constitute only a tiny fraction, / 7 , of the cosmic ray 
flux. When compared with the flux of cosmic rays around the TeV {Ick — 2.582 x 
10 _11 (£'/TeV)~ 2 ' 7 cm" 2 s" 1 sr" 1 MeV -1 [31]) and neglecting gamma attenuation effects 
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one gets the upper limit 

Actually, the attenuation of 7's on the EBL cuts the growth of / 7 , which never exceeds 
10 -5 . Note finally that our extrapolation is consistent with existing observational bounds 
on / 7 at the TeV scale [3l] . 

3. Tracers of the large scale structure 

Since the sources of the CGB are unknown (this is especially true for its high energy 
component of interest here), to predict the anisotropy pattern in the VHE sky we must 
start from some assumption on the distribution of the sources. In the following, we shall 
assume that the sources of the CGB follow the LSS distribution of the matter. Starting 
from Eq. (T5]), the integral flux above the energy E cnt at Earth is simply written as 

F ^^r^ m0w w(E - zh (4) 

where we have defined the window function as 

W(E cut ,z) = / dEg[E(l + z)]e- T ^ z K (5) 

In the limit where the effective cutoff z c « 1, the integrals become almost independent 
on the cosmology, and the above expressions simplify considerably to 

POO 

W(E cut ,z)~ / dEg(E)e-^ E ' z \ (6) 

F(E cut ,h) cc Jdzp a (z,h)W(E cuU z), (7) 

where z is directly proportional to the distance, r ~ cz/Hq. In particular we shall 
assume in the following the power-law g{E) oc E~ 2A , which on the theoretical side 
is consistent with a Fermi-shock acceleration mechanism, and observationally matches 
both the average spectral index of blazar and the spectra of the CGB measured by 
EGRET, see Eq. CO). 



The size of the horizon is roughly 300 Mpc /i -1 at 1 TeV (see Appendix A ) allowing 
safely the use of existing astronomical catalogues like PSCz to predict the gamma signal. 
At lower energy, statistical methods like N-body simulations and the Halo model of 
structures are in principle required to complement existing catalogues, which are not 
deep and/or wide enough. We shall describe these tools in the next section and use 
them to derive the expected statistical properties of the anisotropies at E < 1 TeV. 
These analyses indicate that, at least at large angular scales, most of the non-statistical 
information contained in the surveys can still be retained, allowing one to use the exiting 
catalogues for cross-correlation studies even at E < 1 TeV. 

N-Body Simulations - Our primary tool to compute the statistical properties 
of the matter density distribution is a N-body catalogue. This offers the advantage 
compared to observations of being virtually free of most of the bias affecting astronomical 
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catalogues (no extinction regions, no selection effects, no galaxy morphology and color 
bias, etc). On the other hand, to obtain reliable results one should carefully consider the 
types of bias introduced by the simulation, like the algorithm used, the minimum scale 
resolved, and the physical content and processes treated in the numerical experiment. 

We have performed a set of pure dark matter N-body simulations with a standard 
ACDM universe, a periodic box size of 800 h~ l Mpc, and resolutions of 512 3 and 768 3 
particles using the GADGET2 code [32|,[33]. The initial conditions were computed using 
second order Lagrangian perturbation theory [31] . Placing the observer at an arbitrary 
point in the box, we have computed the density on logaritmically spaced spherical 
shells. A smoothed density field p(z, n) is reconstructed with an adaptive algorithm 
whose smoothing length varies depending on the local number density of dark matter 
particles identical to the one used in the GADGET2 code. After obtaining the window 
functions for different energy cuts, the sky map of the integral gamma ray flux above the 
cut is calculated by integration along radial lines of sight, considering a source density 
proportional to p or to p 2 . We have checked that the resolution used in the simulation 
is sufficient to suppress shot noise from undersampling and other numerical artifacts, by 
comparing the power spectrum extracted at different resolutions. The resulting maps 
F(E CVLt ,n) are obtained on a Healpix grid, facilitating the use of standard tools from 
CMB physics to compute the final power spectrum and multipole coefficients C\ [35J. 
For the final analysis we used 2048 shells, with 12 x 1024 2 pixels on each shell or a 
density grid with almost 25 x 10 9 pixels. This makes it possible to reconstruct the 
power spectrum reliably up to I = 3071. 

Halo Model - In order to check the results of the N-body simulations, we compare 
the Cj's obtained in the case of linear correlations with the predictions using the 
Halo Model [23l [24], [25] . This model is known to provide a fast and efficient semi- 
analytical recipe for describing the clustering of dark matter halos, and of their evolution. 
Comprehensive reviews of the model can be found in [361 137] . Assuming that the 
gamma emitters follow the distribution of the halos, the angular power spectrum of 
the adimensional gamma flux anisotropy [F(E CVLt , h) — (F(E cnt ))] / (F(E cut )) can be 
calculated with the 3D matter power spectrum P(k) and its normalization provided 
by the model. In the flat-sky limit and the Limber approximation, it simply reads 



where the pre-factor N w = J drW(E cnt ,r) represents the contribution from the mean 
intensity (F(E cut )). 

A comparison of the Halo Model and N-body spectra for various energy cuts is 
shown in the top panel of Fig. [2J To properly normalize the flux resulting from the 
N-body catalogue we fix the monopole Cq = 4ir, where the mean flux is derived by 
integrating Eq. (j3j) over the average source density found in the simulation. We see 
that the agreement is pretty good at 100 < I < 1000, while unsatisfactory at very 
low and high /. The former disagreement depends on the breakdown of the flat- 
sky (see e.g. Ref. [38]) and the Limber approximation, while the latter depends on 





(8) 
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Figure 2. In the top panel we show analytical angular spectra (dashed lines) versus 
N-body spectra (solid lines) for energy cuts E cut = 1, 5, 10 TeV (from bottom to top), 
in the hypothesis of linear correlation with matter. In the bottom panel wc show 
N-body spectra for a quadratic correlation. In both cases, the dotted line shows the 
slope of the approximately linear rise of the first multipoles (I <C 100) of the power 
spectrum. The flux from the N-body simulation is normalized by fixing the monopole: 
C = 4tt. 



the inadequacy of the Halo Model to properly take into account non-linearities for 
wavenumbers k > 10 /iMpc" 1 . The biggest disagreement is at 10 TeV, for which the 
window function W is of width ~ 100 h~ x Mpc and the non linear scales are most 
important. In the N-body catalogue we can only account for the flux coming from 
distances up to 350 h^ 1 Mpc. At relatively low energies (E cut < 2 TeV) there are essential 
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contributions to the flux on small scales/high / from distances further away than that, 
while negligible contributions to the flux on large scales/low I. At small angular scales 
and lower energies the simulation also suffers of residual shot noise, seeable in the break 
in the slope at high I, that makes the comparison biased. 

In the lower panel, we present the angular spectrum for the hypothesis of a quadratic 
correlation with density. Both in the linear and the quadratic case, the three most 
interesting features are (i) the power law slope at low scales, which is approximately 
energy independent, as indicated in Fig.[2]by the dotted purple line, (ii) the presence of a 
peak (deriving from the peak of the 3D matter power spectrum), and (hi) the increasing 
relative amplitude of the anisotropics with growing energy cut. At relatively large 
scales, the anisotropics are contributed by the near structures, which are therefore poorly 
affected by the shrinking of the pair-production energy loss horizon. The isotropic or 
small-scale varying component of the flux has instead a significant contribution from far 
objects, which are cutted away at high energies. The main difference between the linear 
and the quadratic case is clearly in the normalization and the slope of the spectrum. 
In the latter scenario the intensity of the anisotropies is enhanced, and the slope is 
steeper. The slope can be estimated even at large scales/low 1, and it is a signature of 
the source correlation with the matter density accessible to observations. For a linear 
density correlation we find /(/ + 1)C/ oc I while for a squared density correlation we find 
l(l + l)Ci oc I 1 - 75 . Also interesting is the prominent shift of the peak to very large I in the 
quadratic case, although it is unlikely that this feature may be observed experimentally. 

We turn now to study how shells at different distances r contribute to the various 
multipoles. Using the predictions of the Halo Model, in Fig. [3] we plot the integrand 
r~ 2 W 2 (E cnt ,r) P(k = l/r) versus r for the two energy cuts E cut =100 GeV and 1 
TeV at different Vs. Note that even for E cnt = 100 GeV, the Q's at / = 10 - 100 
mostly depend on the contribution from structures within a radius of ~ 600 h~ l Mpc. 
Interestingly, the dominant contribution to the first multipoles at I = 1 — 10 comes from 
within a distance of only ~ 300 h^ 1 Mpc. Within such distances astronomical catalogues 
with a large field of view exist (e.g. the 2MASS and PSCz catalogues), providing the 
actual distribution of matter in the universe, and not only a statistical information 
on the density field. Thus, one may aim at the study of the pattern of the first few 
observables a; m 's or, equivalently, of the anisotropy map smoothed to the appropriate 
resolution. Note that, the amplitude of matter fluctuation has a measured maximum at 
a wavenumber (scale) corresponding about to k — 0.02 h Mpc -1 (A = 300 h^ 1 Mpc) [39J. 
Beyond this distance the low I multipoles start to receive decreasing contributions and 
converge rapidly to their asymptotic value. This is a well known expectation, although 
only recently detailed study e.g. of the convergence of the dipole has been performed 
[10] . Actually, as illustrated in Fig. [3] the presence of the window further accelerate the 
convergence of the low / < 10 sky pattern. This point is further illustrated with the 
help of the N-body simulation. In Fig. H] we plot the cumulative contribution to the C/'s 
in function of the maximum distance used in the catalogue, i.e. derived from Eq. (JTj) 
when cutting the integral defining E(E cut ,h) at a distance R. The linear and quadratic 
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Figure 3. Integrand contributions (linear in logr) for the multipolcs I = 
10, 100, 500, 1000 (respectively, black, red, blue, green) and E cut = 100 GeV, 1 TeV. 
Normalization is arbitrary. 

cases are shown, assuming E cut =100 GeV (at higher energies, the effect is even more 
pronounced). Indeed, the simulation confirms that the convergence to the asymptotic 
value is faster at low /'s than at the higher ones. Effectively, for I < 10 — 20 it is 
meaningful to look for correlations between the gamma sky and known local structures 
well below the TeV scale. Of course, the argument can be turned around: by masking 
structures in the nearby sky from the map, the intrinsic low-/ multipoles would be largely 
suppressed, and one may thus access more easily other signatures as the cosmological 
Compton-Getting dipole [28] . 

4. Sky maps 

In the recent years, modern galaxy surveys like the Sloan Digital Sky Survey [4"T1 14"2] 
have greatly improved our knowledge of the distribution of galaxies at large scales, 
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Figure 4. 100 GcV angular spectra when the N-body catalogue is used only until 
a maximum distance R as function of R. Shown are the linear (top) and quadratic 
(bottom) density correlation case. 

revealing a typical foam-like pattern of "filaments and walls" of galaxies around large 
cosmological voids. These very high quality data are not well suited for our analysis 
due to the small fraction of sky surveyed, while we need to perform comparisons with 
the large sky field of view of GLAST and of EAS instruments. In this respect, a fair 
compromise is offered by the IRAS PSCz catalogue [13]. 

The PSCz catalogue contains about 15,000 galaxies and related spectroscopic 
redshifts with a well understood completeness function out to z ~ 0.1. In the limit 
of uniform emission, above E = 100 GeV the majority of the CGB flux is expected 
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Figure 5. Equatorial density sky maps from the PSCz catalogue for the linear (left) 
and quadratic (right) density correlation and for E cut = 100 GeV and 1 TeV. The color 
scale is linear and the average flux outside the mask of the PSCz is normalized to 1 so 
that to represent adimensional maps. The mask of the PSCz survey is enclosed by the 
thick grey contour. 

to come from within this distance. The sky coverage of the catalogue is about 84%; 
the incompleteness is mainly due to the so called zone of avoidance centered on the 
galactic plane and caused by the galactic extinction and to a few, narrow stripes which 
were not observed with enough sensitivity by the IRAS satellite. These regions are 
excluded from our analysis with the use of the binary mask available with the PSCz 
catalogue itself. However, the mask region does not represent a major limitation for 
present analysis, since the galactic emission at low galactic latitudes outshines the CGB 
anyway, prohibiting any CGB analysis in this region. 

We closely follow Ref. [16] for the treatment of the selection effects of the catalogue, 
parameterized via the selection function <p{z). Differently from Ref. [16], we do not sum 
source by source to obtain the final gamma sky map; instead, we use the approach of 
constructing a smoothed density field p(z, n) through the same adaptive algorithm we 
employed in the analysis of the N-body simulation. This method has the advantage of 
efficiently suppressing the intrinsic catalogue shot noise at high redshifts, thus allowing 
us to use the catalogue sources until z = 0.1 ~ 300 Mpc h~ x . We use Eq. PJ, while 
replacing W(E cut , z) — > W(E cut , z)/</>(z), to take into account selection effects. To study 
the sensitivity to the bias of the gamma sources with respect to the baryonic density, 
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Table 1. The characteristics of the experiments considered in our estimates. The 
fraction of the sky observable by a given experiment / s k y is needed for our estimate of 
the errors, see |Appendix B[ 



we shall consider both values a = 1 and a = 2, i.e. a linear and quadratic correlations. 
The latter is in particular expected for dark- matter annihilation models [H]. 

In Fig. [5] we plot the resulting maps from the PSCz catalogue in equatorial 
coordinates for the linear and quadratic density correlation and for E cut = 100 GeV 
and 1 TeV. Let's discuss first the linear case. For the case of the map with E cut = 100 
GeV, modulo the "hole" due to the mask the pattern is quite isotropic with some hot 
spots like from the Virgo and Perseus Clusters. Other structures which appear are 
the Shapley concentration and the Columba cluster (for a key of the local cosmological 
structures see [16] or [15]). From the color scale can be seen that the anisotropy could be 
also a factor of 2 in the hot spots, while is generally of the order 10% at larger angular 
scales as we will show in the following with a multipole decomposition analysis. Given 
the limited statistics of GLAST at high energies, the TeV map is of interest especially 
for the EAS gamma detectors like MILAGRO. We see in this case that the nearest 
structures, forming the Super-galactic Plane, dominate. Of course, from the Northern 
emisphere (where all the present or planned EAS instruments are located) only the 
upper part of the map is visible. Here, the Virgo Cluster and the Perseus cluster offer 
the strongest anisotropy. 

For the quadratic case the change is quite evident: the effect of the quadratic 
correlation is to give more power to the nearest structures (the Virgo and Centaurus 
cluster) that, in fact, are almost dominating the map with fluctuations exceeding even 
10 times the average and appearing almost like point sources. It is instructive to look 
at the case of Shapley concentration at z ~ 0.04 that gives an important contribution 
in all the linear cases but almost disappears in the quadratic maps. The anisotropy in 
the quadratic case is then much more pronounced than the linear case and should be 
easily detectable. 

5. The potential of forthcoming instruments 

In order to estimate the chances of detection of the structures previously described by the 
current or next generation of instruments, one may proceed as follows. First, given the 
specifics of an experiment (in particular its field of view, effective area and background 
rejection capability) one calculates the expected number of events and of misidentified 
cosmic ray background events (iV 7 , N C r) under different assumptions for the E cut and 
the EBL. The events falling in the mask region must be subtracted from the (iV 7 , A*cr) 
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data used for the analysis, so that the incomplete sky coverage is taken into account. One 
then generates M samples of iV 7 events from the VHE 7 map and A^cr from an isotropic 
one. These M realizations may then be used to perform quantitative statements, like for 
example the confidence level with which a given instrument is expected to distinguish 
between the structured sky we predict and an isotropic one, or between linear and 
quadratic correlation scenarios, or between scenarios characterized by a different model 
of CIB. This would be essentially the generalization of the approach followed in [16J. 
Here, however, for the sake of clarity we prefer to develop a simplified analysis in terms 
of the multipole coefficients, ai m , with an analytic estimate of the expected errors due 
to shot noise effects (see Appendix B ). Performing full Monte Carlo analyses for a few 
relevant cases we checked that the errors thus estimated basically agree with the ones 
correctly calculated. In this approach, first we fit the sky maps obtained outside the 
mask via the harmonic expansion 

'max — 10 

F w (E cnt ,h)= a lmYi m (n) ■ (9) 

lm 

As a second step, by knowing the statistics and the characteristics of the detector one 
can estimate the expected errors on the ai m and eventually on the derived coefficients 



Cz's (see Appendix B for details). 



We limit the fit to I < Z max = 10 for several reasons: (i) we have shown that the 
large scales are mainly sensitive to the local universe, for which our predictions are 
robust and deterministic; (ii) the lower the /, the higher the signal to noise ratio is, 
which increases the chances of detection (see below); (iii) the incomplete sky coverage 
due to the mask is affecting our results, and should be taken into account; however, as 
long as we restrict the analysis to the large structures (compared to the size of the mask 
cut) the bias is small. Quantitatively, the method is robust against variations in / max 
as long as the related angular scale 9 min = ir/l max is greater than the typical angular 
extension of the mask cut. More importantly, the purpose of the harmonic analysis is 
to reliably assess the sensitivity of a given experiment; to this aim a small bias in the 
a; m 's is acceptable (we do not want to estimate the true a; m 's, but the detectability 
of them, even of the biased ones). Of course, in the analysis of the real data it will 
be preferable to directly look for cross correlations with the full maps like the ones in 
Fig. H2I whose electronic version is available from the authors upon request. Note that a 
naive a cross-correlation between the gamma sky and the LSS catalogue would not take 
into account the relevant window function effects we have highlighted in this work. 

In the following, we shall consider two kinds of instruments: Satellite-based 
missions, in particular GLAST [16] (to be launched in the Fall 2007), and extensive 
air shower experiments like TIBET |47] . ARGO [IE] (the assembly of which started 
in October 2000 and is presently being completed), or MILAGRO [19]. The latter 
experiment is taking science data since 2001, and since 2004 in its outrigger hardware 
upgrade. In the summer 2007, a transition of MILAGRO to the miniHAWC array 
is expected to start; this is the first step towards the next generation EAS planned 
observatories like HAWC [50]. Unfortunately, ACT instruments like HESS or MAGIC 
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Figure 6. The coefficients ai m up to Z max = 10 calculated from the PSCz gamma 
maps of Fig. [5J The shaded band shows the l-cr shot noise error given by Eq. (|B.6p : in 
the right panel the inner shaded region refers to HAWC, the outer one to MILAGRO. 
We report the predictions for both the linear and quadratic cases. 



are not well suited for such kind of searches, given the small field of view, the low 
duty cycle, and the relatively high impact of the variability of the instrumental and 
atmospheric conditions on the rate of diffuse signals. 

Denoting by I*f(E) the extrapolated EGRET flux which takes into account 
attenuations (see Section [2]), one can estimate the number of events, iV 7 , above the 
energy E y to be collected during the time t as 

POO 

iV 7 = t ■ g cut ■ DC ■ n iov -f m dE A eS {E)I 7 {E) , (10) 

where: DC is the duty-cycle of the instrument; Q{ ov is the solid angle of the field of 
view; / m < 1 is the useful fraction of the sky due to the presence of the galactic mask; 
<7 cu t is the fraction of 7's passing the actual cuts; A e g(E) is the effective collecting area 
of the instrument (averaged over the field of view of the instrument). In the following, 
we shall assume f m = 0.84 due to the mask in the PSCz catalogue, but note that future 
redshift catalogues like 2MRS will have higher / m . Eq. ( iTUl) assumes a quasi-isotropic 
7 sky, which may be violated to some extent at the multi-TeV energies of interest for 
EAS detectors. Even in this case, right ascension anisotropics would not affect the 
estimate, and only large latitude anisotropies might affect A^ 7 by a factor of 0(1). This 
is acceptable enough since we shall only perform a parametric study of the performances 
of an EAS observatory. Analogously, the CR background can be estimated as 

N C R = t- h cut ■ DC ■ Q iov -f m dE A cS (E)I C r(E), (11) 

where now h cut is the fraction of hadrons passing the cuts. Note that we consider the 
same area for CRs as for 7's, although a differential performance of the instrument may 
be taken into account by properly rescaling the factor h cut . The typical parameters we 
shall use are taken from existing literature, and reported in Tab. [TJ Note that GLAST 
is expected to have an excellent background identification, so that only cosmic rays 
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Figure 7. The coefficients Ci up to l ma . x = 10 calculated from the PSCz gamma maps 
of Fig. We shown the level of shot noise [see Eq. (|B.5[) ] expected from four years of 
GLAST and one decade of HAWC. We report the predictions for both the linear and 
quadratic cases. 



in the amount of ~ 6% of the gamma flux pass the cuts. On the other hand, EAS 
experiments have a poor rejection capability (some of them like TIBET jl?] have none), 
which increases typically the gamma content of the diffuse flux by no more than one 
order of magnitude. Therefore one should keep in mind that even after gamma/hadron 
separation, the anisotropics of the gamma sky have to be identified against a quasi- 
isotropic background which is ~ 10 4 larger than the gamma flux. The ultimate limitation 
in detecting anisotropies in the gamma sky with EAS observatories is expected to come 
from the understanding of the intrinsic anisotropy in the CR background. We shall 
come back to this point in the conclusions. 

In Fig. [6] we report the coefficients a; m 's up to Z max = 10 calculated from the 
PSCz gamma maps of Figj5l with the errors estimated according to what reported 
in Appendix B The GLAST mission is expected to last 5+5 years (so we plot a 
realistic 4 years exposure), while EAS instruments have longer run times, so we plot the 
expectations for a decade of collecting time by MILAGRO, or by the proposed project 
HAWC. GLAST should be able to detect some structures above 100 GeV at the 2a level, 
even if the correlation with matter density is only linear. For a quadratic correlation 
one expects a more robust detection, and possibly even hints for anisotropies at higher 
energies^. On the contrary, instruments like MILAGRO may find hints of structures 
(at the 1 a level, see gray band in the right panel of Fig. [6]) only if correlations are 
quadratic or in any case strongly biased with overdensities. As a technical remark, 
note that the performances of MILAGRO above the TeV were estimated by using an 
effective area A^E^) ~ 10 cm 2 (see e.g. [50]). A proper treatment should take 



§ Our estimate does not include the fraction of the CGB measured by EGRET which may be resolved 
by GLAST. If this is removed, our predictions should be rescaled accordingly. 
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into account the energy- dependence of the area, and calculate the expected sky maps 
for MILAGRO weighting accordingly the integral maps. Given the limited chances 
of this instrument to detect the features we have described, we consider this simple 
estimate sufficient to illustrate our point. It is worth to stress that for an EAS 
detector the error on the aim's scales as \/Ncr/N^. Therefore the reduction of the 
shot-noise error goes like (t ■ A c r)~ 1 I 2 (both Ncr and N y grow linearly with t ■ A c g), or 
equivalently as y/h cut / g cut : improving the exposure is equally important as improving 
the gamma/hadron separation capability. A simple inspection of Fig. [6] reveals that for 
a realistic detection of the features in the VHE sky one would need the improvement 
in effective area planned to be reached by instruments like HAWC [50] (see inner green 
band in the right panel of Fig. [6]). An instrument like ARGO [H] is expected to have 
performances in between MILAGRO and HAWC, and may have some chance especially 
if a significant improvement in hadron rejection can be made. Also, note that due to 
their altitude HAWC and ARGO have a significant acceptance of sub-TeV events. While 
the gamma/hadron separation is less efficient at lower energies, the higher statistics may 
help in revealing these structures. 

In Fig. [7] we plot the expectations for GLAST and HAWC in terms of the C\ 
coefficients. At the large scales we are focusing on cosmic variance makes any detection 
of the C{s challenging even when the corresponding easily detectable. This 

proves the importance of the deterministic nature of the expected anisotropics in the 
flux. 

Since our predictions are shaped by the nearby universe, for a fixed background 
the absolute value of the detected anisotropy in principle measures the index of the 
correlation of gamma sources with respect to the matter, as clearly shown by the 
comparison of the top and lower panels in both Fig. [6] and Fig. On the other hand, a 
degeneracy exists with the intensity of the infrared background. While an increase in the 
CIB results in an increased absorption and then in a lower statistics, at the same time 
the horizon for gammas is shrinking and the intensity of the anisotropy increases too so 
that the signal to noise ratio remain almost unchanged. In Fig. [8] we compare the two 
spectra resulting from an increase in the energy cut and an increase in optical depth in 
order to maintain the same collected statistics or, equivalently, the same noise level Cn- 
Indeed, at low energies (E ~ 100 GeV) the change in optical depth only slightly reduces 
the statistics so that the degeneracy is almost perfect. At higher energies (E ~ 1 TeV) 
the change in statistics is more pronounced and the degeneracy is only partial. A suitable 
energy cut that gives the same S/N ratio can still be found, while some differences in 
the multipoles are now visible, although well within typically expected errors. However, 
both the study of single sources and of the energy spectrum of the CGB should pin down 
the remaining uncertainty on r, and the corresponding degeneracy should eventually be 
broken. GLAST, in particular, will put strong constrains on the intensity of CIB from 
the study of 0( 10000) blazars expected to be detected [55] . 

Finally, let's note that these estimates are somewhat conservative: summing the 
power at different Ts may favor the detection (see e.g. [13]), and cross-correlating 
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Figure 8. Upper panel: the coefficients Ci up to i max = 10 calculated from the PSCz 
gamma maps of Fig. [5] as in Fig. \7\ assuming the fiducial value of the CIB optical 
depth and an energy cut E cut = 120 GcV (diamond symbols, black errors bars) and 
a CIB value twice the fiducial one (see Appendix A) and E cut = 100 GcV (triangle 
symbols, red errors bars). The value E cut =120 GeV is chosen so that the noise level 
(dot-dashed line) is the same for the the two CIB values used. Lower panel: same as 
above for E cut = 2 TeV (fiducial CIB value, diamond symbols, black errors bars) and 
E cut — 1 TeV (twice fiducial CIB value, triangle symbols, red errors bars.) 



directly with the maps we have produced would eventually rely on the whole information. 
Indeed, to some extent also the information at small scales can be exploited. In 
particular Fig. [5] shows that the anisotropy expected from small regions of the Virgo 
cluster and along the Super-Galactic plane is significantly above the average flux. 
Therefore, once the data on the CGB become available, reassessing the problem with a 
detailed study will be mandatory. 

6. Summary and Conclusions 

The universe is pervaded by diffuse backgrounds of low-energy photons, of cosmological 
origin like the CMB or due to stellar activity, like the optical and infrared background 
(CIB). This extragalactic background light makes the universe opaque to energetic 
7's. The most energetic part of the gamma-ray background (CGB) is thus of primary 
importance for high energy astroparticle physics, since it acts as a "cosmic calorimeter" . 
Besides telling us about the integrated history of the most powerful astrophysical 
accelerators of the universe, it may constrain non-standard physics taking place at high 
energy, even much higher than the GeV— TeV scale. As an example, we remind that 
the most stringent limits on the decays of superheavy particles coupled at the tree-level 
only to neutrinos come from the observed diffuse extragalactic 7-ray flux |56j . 

In this work we have studied the anisotropy pattern of the CGB in particular in 
the very-high energy regime, beyond about 100 GeV. Due to the onset of the pair- 
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production losses of VHE photons on the CIB, most of the flux is coming from local 
structures, within z < 0.1. Especially at the largest angular scales, the pattern and the 
amplitude of the anisotropies are almost independent of the source energy spectral shape 
and of the cosmological model: modulo the magnitude of the CIB, the key parameter 
in shaping the signature is the degree of the correlation of the gamma-emitters with 
the known matter density field in the nearby universe (i.e. their bias). For example, 
unless substructure emission dominates, dark matter annihilation models for the origin 
of most of the CGB predict a strong (quadratic) correlation of the flux with the matter 
density, which should clearly manifest in the forthcoming observations. It is interesting 
to note the nice complementarity of observations made in the few GeV range with the 
ones in the VHE regime: for example around the GeV the statistics is large and there is 
no uncertainty introduced by the absorption onto the CIB; on the other hand, the VHE 
window is promising since, despite the limited statistics, the anisotropies are larger 
and the predictions are less dependent e.g. on the redshift evolution of the sources. 
Moreover, the "filtering" of the emission of far sources allows one to use more reliably 
the deterministic information given by large scale structure catalogues, going beyond 
merely statistical observables. 

Starting from the PSCz astronomical catalogue, we produced maps of the VHE 
gamma sky and estimated the potential of the satellite mission GLAST or of extensive 
air shower observatories to detect these features. The GLAST mission should be able to 
detect a significant correlation of the diffuse gamma-ray emission with the forecast maps 
presented in this work, providing an important complementary observable to constrain 
the emission models. This is especially true if an exotic contribution from dark matter 
annihilation is relevant. Of course, before claiming an unambiguous detection of dark 
matter in the CGB, detailed particle physics models and a proper foreground analysis are 
required. Indeed, to some extent strongly biased sources may mimick dark matter-like 
features. Some early investigation of this issue are however quite promising [22] . 

EAS experiments are instead limited by the scarce cosmic ray rejection capabilities, 
and only the next generation of instruments like HAWC (or maybe already ARGO) may 
have real chances to achieve the needed sensitivity. For EAS experiments, there is a 
further remark. The Super-Kamiokande experiment has recently detected an anisotropy 
at the level of fewxlO -4 in the cosmic rays around 10 TeV, from a sample of about 
2 x 10 8 muons [52] (MILAGRO has also detected this effect, as mentioned in |53j). 
At similar or higher energies, the TIBET collaboration has reported the detection of 
several anisotropies at the ~ 0.1% level in the cosmic ray flux, probably associated 
with galactic sources and/or galactic transport [54]. While the exposure needed to 
reveal the features we have discussed so far is within the reach of the next generation of 
EAS instruments, even assuming an excellent control over experimental spurious effects, 
the ultimate limitation in detecting these signatures comes from the understanding of 
the intrinsic anisotropy in the CR background. Therefore, an efficient gamma/hadron 
separation is not only necessary to enhance the statistical significance of the point- 
like or diffuse gamma ray sources observed by EAS instruments, but also to control 



The Signature of LSS on the VHE Gamma-Ray Sky 



20 



systematics. In particular, reversing the gamma cut and thus enriching the sample in 
hadronic showers may help identifying and removing non-gamma anisotropies. 
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Appendix A. Gamma Propagation 

Here we report some details on the models used for the extragalactic background light 
(EBL), and on the technique to account for absorption effects in the propagation of 
photons. The main component of the EBL is the cosmic microwave background, the 
spectral number density of which is well known to obey a black-body spectrum 

1 e 2 
7T 2 (he) 3 exp (e/T£ MB ) - 1 

where T<2 MB = 2.73 K=2.35xlO _4 eV is its present temperature. For the 
Infrared/Optical Background (CIB), which is the main source of gamma absorption, 
we use the simple parametrization 

f 5.42xl0 11 eV~ 1 cm^ 3 cxp( ff F) _ 1 A G 200 ^ 2000 /an, 
ncm(e) = < 7.4x lO-VV^cm-V 2 - 295 A G 6.0 -=- 200 /an, (A.2) 

[ rxlO^eV^cm- 3 ^ A < 6.0 /mi, 

where Tp = 13.6 K; and the A — e conversion factor is given by A(/im) = 1.24/e(eV). 
This is a conservative estimate of the CIB consistent with data and constraints reported 
in [29]. 

The time evolution of the backgrounds is obtained simply by redshifting, i.e. 



n(e, z) — (1 + z) 2 n 



€ 



This result is exact for the CMB that is a truly primordial background, while for CIB 
one should in principle perform a simulation of the star formation and dust clustering 
that produced this background in the recent past. However accurate modelling of this 
process suggests that most of the background formed in a burst at z ~ 4 — 5 near the 
peak of star formation rate, so that for our purposes (propagation till z ~ 0.2) simple 
redshifting is quite accurate [29] . 




Figure Al. Critical optical depth r = 1 in function of 7-ray energy and redshift. 
The low CIB case (dashed line) corresponds to Eq. (|A.2|) . while the high CIB case 
corresponds to twice that value. 



In principle, the interaction of high energy photons with the EBL is a complex 
process in which the e~e + created through pair production interact again with 
background photons via inverse Compton scattering producing new high energy 7's. We 
have developed our own code for calculating the evolution of gamma cascades, along 
the lines of Ref. [57] (which we refer to for the details, including the general formalism 
of cascade propagation). Luckily, however, the complicated regime where the cascade 
of high energy e + , e~, and 7 has to be followed in details for calculating the final energy 
distribution is only attained for energies > 10 15 eV. In the TeV range relevant for gamma 
astronomy the development of the shower via secondary particle dynamics can be safely 
neglected, and the observed flux can be calculated by simply considering an energy- 
dependent depletion factor for the spectrum in addition to the effect of redshifting. 
Thus, if the spectrum at the source is of the form g{E), the observed spectral shape of 
the signal will be 

g{E 7 ) oc g[E 7 (l + z)]V{E 7 ,z), (A.4) 

where E 1 is the energy we observe today and V(E^,z) is the probability for a photon 
emitted at redshift z to survive without interacting till now, when it reaches us with 
energy E 1 . This probability is written as 

V(E^z) = e- T{E '" z \ (A.5) 

where the optical depth r is 

r(E y ,z) = f^ ——^—— J den(e, z') J dfi^^a PV (E^l+z'), e, //), (A.6) 
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that is the rate of pair production of photons on the EBL integrated over the time 
during propagating from redshift z to 0. The function <jpp(E, e, /i) is the theoretically 
and experimentally well-known pair production cross section of a photon of energy E 
impinging over a background photon of energy e with a cosinus of the impact angle 
given by \i. In Fig. lAll we show the derived critical t(E 1 , z) — 1 contour that represents 
the gamma redshift horizon as a function of the energy. Note that the EBL is known 
only within a factor of 0(2). As we show in Fig. IA1I with the case of a doubled optical 
depth, there is a degeneracy between the energy horizon and the true value of r. This 
translates into a degeneracy between the anisotropy pattern and the assumed intensity 
of the CIB. However, both the study of single sources and of the energy spectrum 
of the CGB should pin down the remaining uncertainty on this quantity, and the 
corresponding degeneracy should eventually be broken. Given existing uncertainties, 
we find a satisfactory agreement of the function r we computed with more detailed 
studies including the CIB distribution and evolution from simulations of star and galaxy 
formation 1581 . 



Appendix B. Summary of the properties of noise 

In this Appendix we briefly review the noise properties of a discrete poisson process 
on the sphere. The main results can be easily obtained analytically and concern 
the amplitude of the noise variance and the related power spectrum. The map 
resulting from a random realization of A" equally weighted points can be written as 
f(Cl) = 4-ir/NJ2ifi{& — normalized to its mean value, so that to be adimensional; 
then the harmonic expansion coefficients a/ m = / dClfYi m (Jl) follow (apart the constant 
monopole contribution) a gaussian distribution with (ct/ m ) = and 

< m = (4J = % (B.l) 

independent from I, m. From this, the spectrum of the shot noise C\ = £) m \ai m \ 2 / (2/ + 1) 
follows a x 2 distribution with mean 

(Q) = ^ (B.2) 
again independent of /, and variance 

<* = «?>-«*>'=(£)' ^ (B.3) 

These results surely hold in the limit of large statistics; however for our applications 
it is worth testing them also in the limit of very low statistics, say for A" < 1000. 
We performed a set of Monte Carlo simulations to clarify the issue an found that 
the simulations and analytic results are in perfect agreement even in the completely 
unphysical limit of A^ = 10; the analytical results can then be safely used in all the 
cases of interest. 

In general, in addition to the white noise Cj , one has a signal Cf . Moreover, there 
is normally incomplete sky coverage, and additional white noise may be present, as in 
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our case because of the background due to cosmic rays passing the cuts. A generalization 
of Eq. flB.31) then reads 

where / s k y is the fraction of the sky accessible to the experiment (assumed with uniform 
acceptance over this region), and the noise spectrum Cf (including the cosmic ray 
background) is given by 



a 



47r/ sky 



N _ 

7V 7 



^CR 
1 + — ^ 

N 7 



(B.5) 



For the a; m 's, the variance due to the shot-noise plus the background is written as 



°l =Cf = ^# I (l + ^l- (B.6) 

In principle, the complete formulae receive a correction at large I due to the finite angular 
resolution of the experiment that is easily implemented performing everywhere the 
substitution Cf — > Cf exp(Z 2 of /2) cr?, being the angular resolution of the experiment. 
However, since the experimental resolutions are better than a degree, and we limit 
our considerations to the most prominent signatures at large scales, this correction is 
unnecessary for our application. 

We have used these results in estimating the errors reported in Sec. |5j Note that 
what is really measured is always the sum signal+noise: the noise is an unavoidable 
component in any experiment. However, being constant in l,m, the average level of 
the noise can be fitted and subtracted (this is trivial for deterministic predictions since 
( a iL) = 0)- O n the other hand, the error in its determination depends on the sensitivity 
of the experiment and thus the statistics collected, and on the level of background 
rejection. 
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